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The spectral density of various ensembles of sparse symmetric random matrices is analyzed using 
the cavity method. We consider two cases: matrices whose associated graphs are locally tree- 
like, and sparse covariance matrices. We derive a closed set of equations from which the density 
of eigenvalues can be efficiently calculated. Within this approach, the Wigner semicircle law for 
Gaussian matrices and the Marcenko-Pastur law for covariance matrices are recovered easily. Our 
results are compared with numerical diagonalization, finding excellent agreement. 



I. INTRODUCTION 



What started as an approximation to the complex 
Hamiltonian of heavy nuclei has become a very inter- 
esting area of research in its own right. Although the 
statistical properties of random matrices had been tack- 
led before, it was that treatment by Wigner in nuclear 
physics during 1950's which boosted the research of what 
is currently known as Random Matrix Theory (RMT) [l[ . 
The list of applications of this theory has been expand- 
ing ever since, ranging from physics, to computer science 
and finance. Specifically, applications in physics include 
nuclear theory, quantum chaos, localization, theory of 
complex networks, and more (see, for instance, [2| for an 
extensive review). 

From a theoretical and practical viewpoint one of the cen- 
tral quantities of interest in RMT is the spectral density 
of an ensemble of random matrices. While some cases 
have been completely analyzed during the last decades, 
many others have not been fully explored. Consider, 
for instance, the ensemble of symmetric random matri- 
ces whose entries are independently and identically dis- 
tributed Gaussian variables. Among many of its proper- 
ties, it is well-known that its spectral density is given by 
the Wigner semicircle law [3, 13, 3 • Another instance is 
the ensemble of covariance matrices, whose spectral den- 
sity is given by the Marcenko-Pastur law |5|. And the 
list continues. 

Interestingly, the change of introducing sparsity in such 
ensembles, i.e. many entries being zero, complicates the 
mathematical analysis enormously [1, 0, H, M, 03 • Lack- 
ing more powerful mathematical tools, one must rely on 
approximative schemes to the spectral density, e.g. the 
Effective Medium Approximation (EM A), the Single De- 
fect Approximation (SDA) [I,©, [ljj. 
In this work, we tackle the problem of evaluating the 
spectral density of sparse random matrices by using the 
cavity method [ll|, [HJ ■ As we will show, this approach 
may offer new theoretical and practical advantages: from 
a theoretical point of view, it offers an alternative, and 
we believe easier, method to (re)derive the spectral den- 



sity; practically, the resulting cavity equations can be 
interpreted as a belief-propagation algorithm on single 
instances, which can be then easily implemented. The re- 
sulting spectral density is a clear improvement over those 
obtained by approximative schemes. A complementary 
study using the replica method, and emphasizing ensem- 
ble aspects has appeared elsewhere fllij . 
This work is organized as follows: in Sec. [TT]we first men- 
tion how the spectral density can be recast as a problem 
of interacting particles on a graph. The subsequent prob- 
lem is then analyzed by the cavity method in two cases: 
locally-tree like graphs and sparse covariance matrices. 
We derive cavity equations for large single instances and 
check that the dense limit gives the correct results. In 
Sec. IIIII we use the cavity equations as an algorithm to 
calculate the spectral density and compare these results 
with numerical diagonalization. The last section is for 
conclusions. 



II. CAVITY APPROACH TO THE SPECTRAL 
DENSITY 

Consider an ensemble Ai of iVx N symmetric matrices. 
If we denote with {Xf}i=i,...,N as the set of eigenvalues 
of a given matrix A € M. , its spectral density is defined 
as follows 



Qa{X) 



lf>(A-Af). 



(1) 



The spectral density of the ensemble, denoted as p(X) 
results from averaging qa{^) over M.. 
As it was shown by Edwards and Jones [H, the density 
(fT|) can be rewritten as 



Qa{\) = - lim —Trim 
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log Z A (z) 



where 



Z A {z) 



' N 

n 



dx.; 



(2) 



(3) 
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In writing the expression ([3]) , we have been careless with 
the Gaussian integrals, so that as they stand they are 
not generally convergent; we simply follow the prescrip- 
tion as in [lj) EH > and do not worry unnecessarily about 
imaginary factors, so that we can introduce a Gibbs- 
Boltzmann probability distribution of x, viz. 

P A (x) = _!_ e -^(*>*) ( 4 ) 
Za(z) 

with 

1 N 

Ha{x,z) = - ^2 Xi(zl - A)ijXj. (5) 




FIG. 1: Left: Part of a tree-like graph Qa showing the neigh- 
borhood of node i. Right: Upon removal of node i, on the 
resulting cavity graph Q A , the neighboring sites j,k and I 
become uncorrelated. 



In this way, the spectral density qa(^) is recast into a 
statistical mechanics problem of N interacting particles 
x = (#1, . . . , xn) on a graph Qa with effective Hamil- 
tonian ([5]). By Qa we refer to a weighted graph with 
N nodes and edges for each interacting pair with 
weight Aij, when A^ ^ 0. For later use, we introduce 
the following notation: the set of neighbors of a node 
i will be denoted as di; for a given subset of nodes B 
we define xg = (x^ X£, B , ) with £%,..., t\g\ S B and 
with \B\ the number of nodes of B; k% — \di\ denotes the 
number of neighbors of node i, while c = Y2%=i ^ s 
the average connectivity. 

Note that within this approach, Eq. JTJ for the spectral 
density £u(A) can be rewritten as follows: 

1 N 

i— 1 

where (•••)« denotes average over distribution (HJ. 
In previous works [1, 0, H, M EH i the averaged spectral 
densityp(A) is dealt with by using the replica approach, 
or in (lot Il7| using Supersymmetric methods. For gen- 
eral sparse matrices, it was shown in @, H, [ToL llr3| 
that the analysis of the resulting equations from either 
the replica or the Supersymmetric methods was a rather 
daunting task. To push the analysis further, the authors 
in [l0| resorted to a series of approximative schemes, 
originally introduced in 8] . The simpler of such approxi- 
mations, the EMA, assumes that all nodes are equivalent 
and play the same role [1, [l(| • This type of approxi- 
mation works better the larger the average connectivity 
c of the graph. However, for low and moderate values 
of c it fails to provide an accurate description of central 
part and the tails (see, for instance, [9() and of the pres- 
ence of weighted Dirac delta peaks in the spectral density 
[ill, . Other approximations, like the SDA, also fail to 
give an accurate description of the spectrum. 
To improve our understanding of spectral properties of 
sparse matrices, we tackle the problem from a different 
perspective. Instead of considering the averaged spectral 
density p(A), we note, as shown in Eq. ©, that to cal- 
culate gyi(A) we simply need the local marginals Pi(xi) 



from the Gibbs-Boltzmann distribution Pa{x). The cav- 
ity method offers a way to calculate them [24|. To il- 
lustrate this we consider two cases: the ensemble of sym- 
metric locally tree-like sparse matrices, and the ensemble 
of sparse covariance matrices. 

A. Tree-like symmetric matrices 

Let us start by analyzing the spectral density of sparse 
graphs Qa which are tree-like, as the one depicted in 
Fig. [TJ By tree-like we mean short loops are rare. Due 
to the tree-like structure we note that for each node i, 
the joint distribution of its neighborhood P{xgi) is cor- 
related mainly through the node i. If, instead of the 
original graph Qa, we consider a system where the node 
i is removed (see Fig. [TJ, on the resulting cavity graph 
Q^ the joint distribution P^(xgi) factorizes, i.e. 

PW{*Bi)=l[Pi i) (x t ). (7) 

tedi 

This factorisation, which is exact on trees, is called Bethe 
approximation. On the cavity graph, the set of cavity dis- 
tributions {ij (scj)} obeys simple recursive equations, 
viz. 

P^ j \x i ) = ^^- [dx diXj e^ ieBiXj Au*z T| p(i) (xi) 
Z i J eedi\j 

(8) 

for all i = 1, . . . , N and for all j £ di. Once the cav- 
ity distributions are known, the marginal distributions 
Pi(xi) of the original system Qa are given by 

Pi(xi) = e -^- I dx dl e x ^^ A ^ [] P^(x e ) (9) 

* J ii r- a- 

for alii = 1, . . . , JV. While there is in general no a-priori 
reason to expect cavity distributions to be simple, they 
are for the present system: the set (j5J of equations is 

(i) 

clearly self-consistently solved by Gaussian P i s. Hence, 
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upon assuming the cavity distributions to be Gaussian, 
namely, 



P f\x) = 



1 



2ttA 



(10) 



the set of equations ((5j) is transformed into a set of equa- 
tions for the cavity variances Aj(z), viz. 

1 



A«(z) 



(11) 



for all i = 1, . . . , N and for all j S di. Similarly, by Eq. 
© the marginals Pi{xi) are Gaussian with variance Aj 
related to the cavity variances by 



Ai(z) 



1 



(12) 



Eqs. (fTTj) and (fl2|) comprise the final result. For a given 
graph Qai one iterates the cavity Eqs. (fTTj) until conver- 
gence is reached. Once the cavity variances are known, 
the variances A, are given by Eqs. (|12[) . from which the 
spectral density qa{^) is obtained by 



Qa(X) 



lim 

e^0+ TlN 



N 



^Im[A,(z)] 2 



(13) 



Large c limit: The Wigner Semicircle Law 

To assess our approach, note that from the set of equa- 
tions (fTT|) and (fT2"j) we can easily recover the Wigner semi- 
circle law in the large c limit. By this limit we understand 
that the fe, — > c and c — > oo. and assume that the graph 
is already "infinitely" large [25|]. To do this large c limit, 
we take the entries of the matrix A to be Aij — Jtj /^/c, 
with Jij{= Jji) a Gaussian variable with zero mean and 
variance J 2 . From Eqs. (fTTj) and (| 1 2[) . we note that 

for large c, we have that Ap'(z) = A,(z) + ©(c" 1 ) 
Upon defining 



A 



lim 

c — >oo c 



r * — ✓ 



(16) 



we obtain that 

lim £ A 2 £ A« = lim 1 ^ J ^ = j2A ■ ( 17 ) 



Thus, in the large c limit Eq. (fT2"j) yields 

1 



A = 



z - J 2 A ' 

which gives the well-known Wigner semicircle law Q 



(18) 



p(A) = ^v/4J 2 -A 2 

27TJ Z 



(19) 



»=i 



It is worth noting that equations equivalent to those in 
(fTTj) can be derived by the method described in 1201] and 
can be related to self-returning random walks l2l|. 
Notice also that the set of cavity equations must be solved 
for complex z = \ — ie, so that the cavity variances are in 
general complex, and then the limit e — ► + is performed. 
Instead, we perform this limit explicitly in the cavity 
equations. To do so, we separate these equations into 
their real and imaginary parts and then do explicitly the 
limit e — > + by naively assuming that the imaginary part 
is non- vanishing in such a limit (see discussion below). 



Denoting (a- , b. 



[Re(Ap ) ),Im(A^ ; )], we obtain 



a9 ] = 



\-h<?\a) 



di) 



\~hf\a)) +(h<f\b))' 
(x-h^(a)y+(h^(b))'' 



(14) 



with 



(15) 



Our results are exact, as long as the average connectivity 
c of the graphs considered remains finite in the limit N — > 
oo. 



B. Covariance matrices 



Let us consider now matrices A of the type 



-t 

d i 



(20) 



where £ is an N x P matrix with entries To this 
matrix we can associate a bipartite graph Q% with N + P 
nodes, divided into two sets indexed by i = 1,...,N 
and \x — 1,...,P (see Fig. [2]). A pair of nodes 
is connected if &p ^ 0. We refer to these nodes as x- 
nodes and m-nodes, respectively. We will consider the 
bipartite graph to be tree- like, i.e. many of the entries 

_ x kfj, with 



£i M are zero. We also introduce d = (1/P) J2, 



kfj, = i.e. the average connectivity of the m— nodes. 
Clearly, c = ad with a = P/N . In this case it is more 
convenient to apply the cavity method on the bipartite 
graph To do so we write the effective Hamiltonian 
([5]) as follows 



1 N 1 P 



i— 1 M— 1 

where we have defined the overlaps 



(21) 



(22) 
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FIG. 2: Left: Graph Qa for covariance matrices. Right: Bi- 
partite graph Q$ for the matrix £. For sake of clarity, self- 
interactions in the graph Qa are not drawn. 



equations for the cavity variances 



■ _ I V f 2 i 



(26) 



eedu\i 



Similarly, if we denote with A, the variance of the 
marginal Pi(xi) on the original graph we obtain 



A;(z) 



_ I V £ 2 — 



(27) 



Note that, due to our choice for £ and the relation (|2"0|) 
between the matrices A and the corresponding graph 
Qa (see Fig. [2]) is locally clique-like with self-interactions. 
We work our cavity equations in the bipartite graph Q^. 
Here, the variables x are on the a;-nodes while the vari- 
ables m — (mi, . . . , m p ) are on the m-nodes. Since we 
have two types of nodes, we apply the cavity method 
twice: around a;-nodes and around m-nodes. We define 
Qv (m v ) as the cavity distribution of m„ in the absence 
of a node i, and P^ p \xi) is the cavity distribution of xi 
in the absence of node fi. We find the following set of 
equations for these cavity distributions 



g 2 ZX i 



dm 



(23) 



for all i = 1, . . . , N and [i £ di. Also 



K) 



dx du \i S | ;/;,, 



1 



x n 

iedu\i 



(24) 



for all f = l,...,P and i G dv. Obviously, for the 
marginal distributions Pi{xi) we obtain 



Pi(Xi) 



e 2 * 



drriQi e 



n 



(25) 



for all i = 1, . . . , N. As before, we see from the set of 
equations |23|) and (|24|) that the Gaussian measure is a 

fixed point. Thus, by taking P^{xi) and Q^\m^) to 
be Gaussian distributions with zero mean and variances 
A^ and , respectively, we obtain the following set of 



Large c limit: The Marcenko-Pastur law 

For the sake of simplicity we take the non-zero entries 
of the matrix £ to have values ±1, so that £ 2 M = 1 in Eqs. 
(|26|) and (l27|) . Let us consider the spectral density in the 
large c limit of the bipartite graph . By this limit we 
mean — > d, fej — > c and d, c — > oo while a remains 
finite. As before, the difference between cavity variances 
and variances is 0(c~ r ). If we define 



1 1 d 

A = lim - A e = lim -S^ A t 

d-too d — ' d^oo d Z — ' 



(28) 



eedu 

from Eq. (|27|) we obtain 
A = 



z — a 



1-A 



(29) 



Upon solving this equation for Im(A) we obtain the 
Marcenko-Pastur law [5J] of dense covariance matrices 



1 



.(A) = — ^T2A(^Ti)T(^i? (30) 

+ C (l - a)S(X) . 

with Co = 1 for a < 1 and Cq = for a > 1. A slightly 
different expression is found by Nakanishi and Takayama 
[22I ] , where the difference comes from not considering the 
diagonal terms. This could also be implemented f airly 
straightforwardly to obtain the spectral density as in [22| • 



III. NUMERICAL RESULTS AND 
COMPARISON 



For general sparse matrices, we solve the cavity equa- 
tions numerically and compare the results with exact nu- 
merical diagonalization. We consider again the two cases 
of locally tree- like and sparse covariance matrices. 



X 



X 



FIG. 3: Spectral density of Poissonian graphs with bimodal 
edge weights and average connectivity c = 3. Red Square 
markers are the results of numerical diagonalization with 
N = 1000, averaged over 1000 samples. Blue circles are the 
results of the cavity approach with N = 1000, averaged over 
1000 samples. The dashed line corresponds to the SDA and 
dotted line is the EMA. The inset shows the tail of the spec- 
tral density. 

A. Tree-like symmetric matrices 

To test our cavity equations, we choose Poissonian 
graphs Ga where each entry Aij of the N x N matrix 
A is drawn independently from 

P(A ij ) = ^(A ij ) + (l-^)6(A ij ) (31) 

with c the average connectivity, and n(x) is the distri- 
bution of non-zero edge weights. For the distribution of 
weights we study two cases: bimodal distribution, i.e. 

7r(^) = \S{A 13 - 1) + ^(Aj + 1) (32) 

and Gaussian distribution with zero mean and variance 
1/c. 

For the purpose of fairly comparing later with exact 
numerical diagonalization, we have analyzed the cavity 
equations for rather small matrices. However, we have 
checked that the convergence of these equations is gen- 
erally fairly fast and we are able to evaluate the spec- 
tral density of very large matrices in reasonable time. In 
both, the bimodal and the Gaussian cases, we generated 
matrices with N = 1000. For each matrix we run our 
cavity Eqs. (jlip until convergence is reached and then 
obtain the spectral density from Eqs. I|12p. The result 
is averaged over 1000 samples. For such sizes we have 
also calculated the spectral density by exact numerical 
diagonalization and averaged over 1000 samples. 

The numerical results from the cavity approach and 
exact numerical diagonalization for the bimodal case is 
plotted in Fig. [3] for average connectivity c = 3. We 



FIG. 4: Comparison of cavity equations for e = (blue circles) 
and e = 0.005 (Continuous red line), for Poissonian graphs 
with bimodal edge weights and average connectivity c = 3 
(N — 1000 and average over 1000 samples). The inset shows 
the Dirac delta structure in the central region. 



have also compared our results with the spectral den- 
sity obtained by EMA and SDA (see [|, H, H| for details 
about the approximations). As we can see, our results 
are a clear improvement over the EMA and SDA results, 
as they are in excellent agreement with numerical diag- 
onalization. Even the tail of the spectrum, usually not 
obtained with these approximation, (see inset of Fig. [3]) 
is well reproduced by our approach. 
It is well known that the spectrum of these type of en- 
sembles contains a dense collection of Dirac delta peaks 

nana, 

which are not fully captured by the previous 
approximations. Without a prior analysis, one wonders 
how the cavity equations can be used to obtain such 
contributions. A practical way out is to reconsider the 
limit e — > 0, by leaving a small value of e in the cavity 
equations, which implies approximating Dirac deltas by 
Lorentzian peaks. In Fig. U]we have rerun the set of eqs 
(jlip and (|12p with a small value of e. The Dirac delta 
contributions, whose exact positions within the spectral 
density is discussed in [TH, are now clearly visible. A 
more detailed study on this issue within the context of 
localisation can be found in [l3j . 

In Fig. [5] we plot the results of both numerical diago- 
nalization and the cavity method when tt(x) is a Gaus- 
sian distribution with zero mean and variance 1/c. Once 
again, our results are in excellent agreement with the nu- 
merical simulations. 



B. Covariance matrices 

We have also analyzed numerically in the case of sparse 
covariance matrices. Here the entries ^ of the N x P 
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FIG. 5: Spectral density of Poissonian graphs with Gaussian 
edge weights and average connectivity c = 4. Red Square 
markers are for the results of numerical diagonalization with 
N — 1000, averaged over 1000 samples. Blue circles are re- 
sults of the cavity approach with N = 1000 and average over 
1000 samples. 



0.25 - 



p(X)0.15 




0.05 



FIG. 6: Spectral density of covariance matrices with N = 
4000, d = 12, a = 0.3. Average over 1000 samples. Red 
Square markers are for the results of numerical diagonaliza- 
tion. Blue circles are results of the cavity approach. The 
dashed line corresponds to the SEMA. 



medium approximation (SEMA), introduced in [lf]l |; here 
we make the same choice of parameters. The inset figure 
shows detail of the tail region of the plot, where the dif- 
ference between the SEMA and the other results can be 
clearly seen. 



IV. CONCLUSIONS 

In this work, we have re-examined the spectral den- 
sity of ensembles of sparse random symmetric matrices. 
By following Edwards and Jones [4j, we have mapped 
the problem into an interacting system of particles on a 
sparse graph, which was then analyzed by the cavity ap- 
proach. Within this framework, we have derived cavity 
equations on single instances and used them to calculate 
the spectral densities of sparse symmetric matrices. Our 
results are in good agreement with numerical diagonal- 
ization and are a clear improvement to previous works 
based on approximative schemes. We have also shown 
that, to account for the Dirac delta contribution to the 
spectrum, one may approximate Dirac delta peaks by 
Lorentzians, by leaving a small value for e [HI]. 
It is well known that cavity and replica methods are 
equivalent (see for instance [l2| for diluted spin glasses), 
so one may wonder in which aspects our work differs 
from the ones presented in 0, H, OH ■ Generally, for 
interacting diluted systems with continuous dynamical 
variables one expects an infinite number of cavity fields 
to parametrize the cavity distributions. The authors in 
0, H, H, [H| decided to tackle such a daunting task by 
resorting to approximations. In this work, we simply 
realize that, for the problem at hand, the cavity distri- 
butions are Gaussian, so that the problem can be solved 
exactly by self-consitcntly determining the variances of 
these distributions. 

In future studies we expect to extend the method pre- 
sented here to the analysis of more general aspects of 
random matrices. 
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